function runfamamacbethreturns_quarterlyreturns_boom_bust(merged,programFiles,Outfile,qtrdummy,ff_fact,boomdummy)

%prepare the variable for output file
outputFormat.Intercept = NaN*ones(500,1);
outputFormat.beta_lag_1q = NaN*ones(500,1);
outputFormat.fill_size_lag = NaN*ones(500,1);
outputFormat.fill_bm_lag = NaN*ones(500,1);
outputFormat.fill_lev_lag = NaN*ones(500,1);

outputFormat.fill_hhi_lag = NaN*ones(500,1);
outputFormat.fill_ol_lag = NaN*ones(500,1);
outputFormat.hhi_ol = NaN*ones(500,1);

outputFormat.fill_op_pft_1_lag = NaN*ones(500,1);
outputFormat.hhi_op_pft_1_lag = NaN*ones(500,1);


%prepare interactions
merged.hhi_ol = merged.fill_hhi_lag.*merged.fill_ol_lag;
merged.hhi_op_pft_1_lag = merged.fill_hhi_lag.*merged.fill_op_pft_1_lag;
outputFormat = struct2table(outputFormat);




%prepare all specifications
Specifications{1} = {'beta_lag_1q','fill_size_lag','fill_bm_lag','fill_lev_lag','fill_hhi_lag'};
Specifications{2} = {'beta_lag_1q','fill_size_lag','fill_bm_lag','fill_lev_lag','fill_hhi_lag','fill_ol_lag','hhi_ol'};
Specifications{3} = {'beta_lag_1q','fill_size_lag','fill_bm_lag','fill_lev_lag','fill_hhi_lag','fill_op_pft_1_lag','hhi_op_pft_1_lag'};


%rearrange data
merged = sortrows(merged,{'GVKEY','datadate'});





%returns are already in percentages i.e., 8% is 8 in the file
if qtrdummy == 1
    DepVar = 'qtrret'
    disp('Note: Running the regression every quarter')
    ff_fact.Mkt = 1+ff_fact.Mkt/1e2;
    ff_fact.Mkt = [nan;nan;ff_fact.Mkt(3:end).*ff_fact.Mkt(2:end-1).*ff_fact.Mkt(1:end-2)];
    ff_fact.Mkt = (ff_fact.Mkt - 1)*1e2;

    loc = find(ff_fact.Mkt > 10);
    if boomdummy == 1
        ff_fact = ff_fact.yyyymm(loc);
    else
        ff_fact = setdiff(ff_fact.yyyymm,ff_fact.yyyymm(loc));
    end
end



forregression(Specifications,DepVar,[],merged,outputFormat,qtrdummy,programFiles,Outfile,ff_fact);





function forregression(Specifications,DepVar,nlag,merged,results,qtrdummy,programFiles,Outfile,ff_fact)
if qtrdummy == 1
    nlag = 2;

    disp('Note: Running the regression every quarter')
    disp(strcat('Lags: ', num2str(nlag)))
    loc = merged.qtrdummy == 1;
    merged = merged(loc,:);

end
counter = 1;
counter1 = 1;
results1 = results;
resultsMean1 = results;
resultsMean = results;
for jm = 1:length(Specifications)
    disp(jm)


    [Coeff, Mn, forMnpool] = performregressions(table2array(merged(:,DepVar)),merged,Specifications{jm},nlag,merged.datadate,ff_fact);
    [ir,ic] = find(isnan(Coeff) == 1);
    Coeff(ir,:) = [];
    Dates = Coeff(:,1);
    Coeff = Coeff(:,2:end);

    tmpdates = floor(Dates/1e2);
    [~,ia,ib] = intersect(tmpdates,ff_fact);
    BoomDummy = zeros(size(Dates));
    BoomDummy(ia) = 1;
    NoBoom = ones(size(Dates)) - BoomDummy;
    % Dates = Dates(ia);
    
    clear CoeffAll temptStat tempSE CoeffAllno temptStatno tempSEno


    %compute the newey west t-statistics
    for j = 1:size(Coeff,2)
        y1 = [ones(size(Coeff(:,j))), Coeff(:,j)];
        [ir,~] = find(isnan(y1) == 1);
        y1(ir,:) = [];
        tempBoomDummy = BoomDummy;
        tempNoBoomDummy = NoBoom;
        tempBoomDummy(ir) = [];
        tempNoBoomDummy(ir) = [];
        [a,~,c] = regress(y1(:,2),[tempBoomDummy, tempNoBoomDummy]);
        nwse = neweywest(c,[tempBoomDummy, tempNoBoomDummy],nlag,0);
        
        CoeffAll(1,j) = a(1);
        temptStat(1,j) = a(1)./nwse(1);
        tempSE(1,j) = nwse(1);

        CoeffAllno(1,j) = a(2);
        temptStatno(1,j) = a(2)./nwse(2);
        tempSEno(1,j) = nwse(2);

        

    end

    tStatAll_1 = temptStat; clear temptStat;
    % CoeffAll = nanmean(Coeff);
    Coeff = CoeffAll; Signif = tStatAll_1; SE = tempSE;
    


    % clear CoeffAll tStatAll;
    Final = [];
    for jj = 1:size(Coeff,1)

        temp = [Coeff(jj,:); Signif(jj,:); SE(jj,:)];
        Final = [Final; temp];

    end

    tStatAll_1 = temptStatno; clear temptStat;
    % CoeffAll = nanmean(Coeff);
    Coeff = CoeffAllno; Signif = tStatAll_1; SE = tempSEno;

        Final1 = [];
    for jj = 1:size(Coeff,1)

        %     temp = [Coeff((jj,:)); Signif((jj,:)); Signif((jj,:))*NaN];
        temp = [Coeff(jj,:); Signif(jj,:); SE(jj,:)];
        Final1 = [Final1; temp];

    end
    
    [results,counter,resultsMean] = arrangeresults(counter,results,Final,Specifications,jm,Mn(BoomDummy==1,:),resultsMean);
    [results1,counter1,resultsMean1] = arrangeresults(counter1,results1,Final1,Specifications,jm,Mn(NoBoom==1,:),resultsMean1);



    


end

%write the results in excel file
check = exist(strcat(programFiles,'\output\'),'dir');
if check == 0
    mkdir(strcat(programFiles,'\output\'))
end
cd(strcat(programFiles,'\output\'))
writetable(results,Outfile,'Sheet',strcat(DepVar,'_Boom'));
writetable(resultsMean,Outfile,'Sheet',strcat(DepVar,'_Boom'),'Range',strcat('A',num2str(3*(length(Specifications)+1))));

writetable(results1,Outfile,'Sheet',strcat(DepVar,'_NoBoom'));
writetable(resultsMean1,Outfile,'Sheet',strcat(DepVar,'_NoBoom'),'Range',strcat('A',num2str(3*(length(Specifications)+1))));


cd(programFiles)


function [results, counter,resultsMean] = arrangeresults(counter,results,Final,Specifications,j,Mn,resultsMean)
    [~,~,ir] = intersect({'Intercept'},results.Properties.VariableNames);
    out = Final;
    St = counter;
    En = counter + 2;
    results(St:En,ir) = array2table(out(:,1));
    
    for k = 1:length(Specifications{j})
        [~,~,ir] = intersect(Specifications{j}{k},results.Properties.VariableNames);
        %         if j == 14
        %             pause
        %         end
        %         Specifications{j}{k}
        %         array2table(out(:,k+1))
        %         disp(ir)
        %          results(St:En,ir)
        results(St:En,ir) = array2table(out(:,k+1));
        resultsMean(St:St,ir) = array2table(mean(Mn(:,k)));
    end
    counter = En + 1;


function [Coeff, Mn, forMnpool] = performregressions(y,xTable,Spec,nlag,dates,ff_fact)
%Inputs
%y: the dependent variable
%xTable: contains all the independent variables
%Spec: the specification for which we want to run the regression
%nlag: number of lags for newey-west

x = [];
for j = 1:length(Spec)
    disp(Spec{j})
    [~,~,ic]= intersect(Spec{j},xTable.Properties.VariableNames);
    disp(xTable.Properties.VariableNames{ic})
    x = [x, table2array(xTable(:,ic))];
end

x = [dates,y,ones(size(y)),x];
[ir,~] = find(isnan(x) == 1);
x(ir,:) = [];

dates = unique(dates);



forMnpool = [];
for j = 1:length(dates)

    loc = find(x(:,1) == dates(j));
    yy = x(loc,2);
    xx = x(loc,3:end);

    tx = [yy,xx];
    [ir,~] = find(isnan(tx) == 1);
    tx(ir,:) = [];
    [ir,~]  = find(isinf(tx) == 1);
    tx(ir,:) =[];


    yy = tx(:,1);
    xx = tx(:,2:end);

    R = rank(xx);

    Mn(j,:) = mean(tx(:,3:end));
    forMnpool = [forMnpool; tx(:,3:end)];

    if R == size(xx,2) && size(tx,1) > 7
        %         disp(size(yy))
        %         disp(size(xx))

        [b,~,~,~,tR2] = regress(yy,xx);
        R2(j,:) = tR2(1);
        Coeff(j,:) = [dates(j), b'];
    else
        if j == 1
            Coeff(j,:) = [dates(j), ones(1,size(xx,2))*NaN];
        else
            Coeff(j,:) = [dates(j) Coeff(j-1,2:end)*NaN];
        end
    end




end


